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Biofilms are antibiotic-resistant bacterial aggregates that grow on moist surfaces and can trigger 
hospital-acquired infections. They provide a classical example in biology where the dynamics of 
cellular communities may be observed and studied. Gene expression regulates cell division and 
differentiation, which affect the biohlm architecture. Mechanical and chemical processes shape the 
resulting structure. We gain insight into the interplay between cellular and mechanical processes 
during biohlm development on air-agar interfaces by means of a hybrid model. Cellular behavior is 
governed by stochastic rules informed by a cascade of concentration fields for nutrients, waste and 
autoinducers. Cellular differentiation and death alter the structure and the mechanical properties 
of the biohlm, which is deformed according to Foppl-Von Karman equations informed by cellular 
processes and the interaction with the substratum. Stiffness gradients due to growth and swelling 
produce wrinkle branching. We are able to reproduce wrinkled structures often formed by biofilms 
on air-agar interfaces, as well as spatial distributions of differentiated cells commonly observed with 
B. subtilis. 

PACS numbers: 87.10.Mn, 87.18.Hf, 87.18.Fx 


I. INTRODUCTION 

From bacterial communities to multicellular tissues, 
three dimensional biological structures emerge through 
poorly understood interactions between mechanical 
forces and cellular processes. Being apparently less com¬ 
plex than multicellular organisms, bacterial biofilms may 
provide model systems for analyzing the interplay be¬ 
tween cellular and mechanical features of three dimen¬ 
sional self-organization during growth. 

Biofilms are bacterial aggregates that form on moist 
surfaces and are sheltered from external aggressions 
by a self-produced extracellular matrix (ECM) made 
of exopolymeric substances (EPS) p], 2 . This matrix 
makes them uncommonly resistant to antibiotics, disin¬ 
fectants, flows and other chemical or mechanical agents 
[3]. Biofilms are involved in chronic infections associ¬ 
ated to biomedical implants HE] and in many industrial 
problems, such as biocorrosion, biofouling, efficiency re¬ 
duction in heat exchange systems and food poisoning 0 - 
:8 . They are beneficial in bioremediation and biocontrol 
0EEQ] or in wastewater treatments m- 

Biofilm structure depends among other variables on 
the bacterial strain and the nutrient source, but is also 
influenced by environmental conditions. Biofilms grown 
in flows mm differ noticeably from biofilms expand¬ 
ing on air-solid or air-liquid interfaces not exposed to 
a shear force must These vary depending on the sur- 
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face, as it is the case for Bacillus subtilis forming pellicles 
on air-fluid interfaces [16] or wrinkled films on agar m- 
Pattern formation is crucial to the development of biolog¬ 
ical systems. In cellular organisms, growth is a complex 
process implicating biochemical and physical mechanisms 
occurring at a variety of length and timescales [18 . On 
one side, genetic programs govern cellular processes, such 
as growth and differentiation. On the other, mechan¬ 
ical properties and physical forces induce macroscopic 
movements of cell populations and chemical compounds 
mum. Detailed models trying to describe biofilm evolu¬ 
tion on surfaces should take into account the interaction 
of these processes. 

Many models have been proposed for biofilm expan¬ 
sion in flows H3I2D, focussing typically on growth due 
to nutrient consumption. The standard geometry con¬ 
sists of a biofilm attached to a solid substratum sub¬ 
merged in a fluid. Growth is limited by diffusion of 
oxygen and nutrients from the surrounding flow into the 
biofilm. EPS production favors vertical spread, pushing 
cells into the flow to improve access to nutrients and oxy¬ 
gen [lj [22] . Depending on the nutrient and oxygen avail¬ 
ability, the bacterial strain and the hydrodynamic condi¬ 
tions, biofilms evolve to form a variety of patterns: mush¬ 
rooms, streamers, mounds, ripples and patchy shapes are 
observed P1|T31P31|22]. The volume fraction of EPS in 
these biofilms is large m • Some models view them as 
cells embedded in a fluid EPS matrix [25 . Experimen¬ 
tal observations and measurements suggest that biofilms 
behave as viscoelastic materials, i.e., exhibit elastic solid- 
like response to short time scale stimuli and viscous fluid- 
like response to long time scale stimuli [26 . Two-phase 
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flow models of biofilm spread in fluids m coexist with 
models that consider the growing biofilm an hyperelas¬ 
tic material [28] . The interplay of biofilm filaments with 
the flow, however, is described as interaction between an 
elastic structure and a fluid H3H3- These models rep¬ 
resent the microbial community as a continuum. To in¬ 
vestigate the influence of interactions at bacterium level, 
‘agent’ based models have been developed. Bacteria are 
described as individuals whose behavior responds to ex¬ 
ternal continuum fields. Different modeling frameworks 
characterize individual bacteria in different ways. Cellu¬ 
lar automata [3T1433] . cellular Potts [34] and individual- 
based [22] [25] |35] representations have been proposed for 
biofilms in flows. We will comment on the general ap¬ 
proaches in Section II. Basic mechanisms incorporated 
in them include biofilm expansion due to cell division, 
phenomenological descriptions of EPS production and 
biofilm erosion caused by the external flow. 

Biofilms grown in other environments show a differ¬ 
ent structure and undergo different processes. Under¬ 
standing their evolution requires different models, that 
must take into account the specific behavior of the se¬ 
lected bacterial strain. We consider here the develop¬ 
ment of B. subtilis biofilms on air-agar interfaces, for 
which detailed experimental evidence is becoming avail¬ 
able damans gsi bzi. Recent observations reveal the 
role of cellular behavior and inner mechanical processes 
on the biofilm constitution. These biofilms extract nu¬ 
trients and water from the agar substratum itself. Cells 
are densely packed, glued together by a small volume 
fraction of extracellular substances [36j- The specific 
autoinducer sequence triggering EPS secretion has been 
identified 0 . EPS production favors water absorption 
and horizontal spread [36] . Localized cell death leads to 
mechanical stress relief through buckling d3. Instead 
of networks of vertical mushrooms or filaments undulat¬ 
ing with the external flow, wrinkled films are observed 
ma na eu. Partial aspects of this picture have been 
modeled. Ref. [36] analyzes accelerated horizontal ex¬ 
pansion due to EPS production and swelling by means 
of thin film approximations calibrated to experiments. 
A reaction-diffusion model for the density of living cells 
combined with equations for the accumulation of EPS 
and waste provides understanding of the cell death pat¬ 
terns in Ref. m- We introduce here a hybrid framework 
that might be suitable to simulate the interplay between 
the mechanical and cellular processes identified so far. 
We couple a stochastic treatment of individual cellular 
activities with a continuum description of elastic defor¬ 
mations. To facilitate the coupling of stochastic cell di¬ 
vision, death, differentiation and swelling processes with 
macroscopic elastic deformations and the description of 
the complex interaction with the substratum, we have 
represented bacteria as cellular automata. Cells respond 
to a number of concentration fields (nutrients, waste, 
secreted substances), and are affected by water migra¬ 
tion processes and elastic deformations resulting of inner 
stresses due to growth, swelling and the interaction with 


the substratum. Our numerical simulations suggest that 
the wrinkle branching observed in experiments may be 
the result of elastic deformations triggered by stiffness 
gradients. 



FIG. 1. (Color online) (a) B. subtilis biofilm grown on a 1.5% 
agar surface in a biofilm-inducing medium. Photograph taken 
after four days. Reprinted with the permission of Cambridge 
University Press from Chai L, Vlamakis H, Kolter R, MRS 
Bulletin, 36: 374-379 (2011), [TJ. (b) In silico film showing 
a corona of radial wrinkles emerging from a central core. 

Typical wild type B. subtilis biofilms on agar show a 
central core with highly connected wrinkles surrounded 
by a corona with radial branches, as in Figure [l] Succes¬ 
sive branching is also reported, see Figure [2] These struc¬ 
tures differ from those commonly observed in biofilm pel¬ 
licles formed on air-liquid interfaces m • Such pellicles 
spread on a space confined by the walls of the vessel con¬ 
taining the liquid. Ref. m assumes an expression for the 
compression stresses and performs a stability analysis of 
the equilibria of a plate model, concluding that a buckling 
instability should appear above a threshold. The pres¬ 
ence of walls stresses the pellicles as they grow, produc¬ 
ing wrinkles. In contrast, the biofilms studied here spread 
on an air-agar interface with no boundaries around. The 
mechanism for wrinkle formation is bound to be differ- 
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FIG. 2. (Color online) (a) B. subtilis biofilm grown on the 
air-solid surface of agar gel containing water and nutrients. 
Reprinted with the permission of Cambridge University Press 
from Wilking JN, Angelini TE, Seminara A, Brenner MP, 
Weitz DA, MRS Bulletin, 36: 385-391 (2011) [H]. (b) In 
silico film showing successive wrinkle branching. 


ent, even if elastic deformations are involved too. Ref. 
[38] uses a neo-Hookean energy to study the elastic de¬ 
formations of a growing circular film moderately adhered 
to a rigid substrate and subject to anisotropic growth. 
Tracking bifurcations in families of analytical solutions 
as the anisotropy degree increases, they predict the ap¬ 
pearance of contour undulations, delamination, buckling 
on the scale of the sample after delamination and, finally, 
a periodic distribution of folds. As Ref. [38 points out, 
such buckling patterns are not experimentally observed 
in B. subtilis biofilms. However, the analysis may be use¬ 
ful to understand delamination and folding phenomena. 
Those models do not consider the microscopic processes 
occurring in a biofilm. The complex interaction with the 
substratum is reduced to attachment using a collection 


of springs. On the other hand, neither energy nor linear 
stability approaches answer the question of how different 
patterns emerge and evolve from one to another. Here 
we aim to compute the development of wrinkled struc¬ 
tures as the biofilm expands horizontally, starting from 
an almost fiat initial state. We will solve time depen¬ 
dent equations coupling the film to the substratum. The 
residual stresses will be computed numerically averaging 
information from growth and water absorption processes. 
This allows to observe how wrinkles develop and how 
they respond to localized variations in the microstruc¬ 
ture. We recall below some experimental observations to 
bear in mind in the sequel. 

During biofilm expansion, B. subtilis populations dif¬ 
ferentiate in several types in response to local environ¬ 
ments created by growth, waste production, nutrient con¬ 
sumption and cell-cell communication [14] . Cells produce 
a signaling molecule, ComX. When a threshold concen¬ 
tration is reached, a percentage of the population of cells 
expresses some genes that induce a change in their pheno¬ 
type, becoming surfactin producers. Surfactin may signal 
other cells (that do not make surfactin) to produce ex¬ 
tracellular matrix. As the biofilm thickens, starving cells 
in the top become inert spores. 

Initially, part of the cells are motile and express flag¬ 
ella allowing to swarm on the surface. As the production 
of ECM increases, most cells lose their individual motil¬ 
ity m- Flagella are downregulated and biofilm spread 
is influenced by ECM production [36]. In early stages 
of the biofilm evolution, dividing cells interact with their 
neighbors and push them. The biofilm becomes thicker 
and the edges spread slowly on the surface. Later in¬ 
crease in the exopolysaccharide concentration raises the 
osmotic pressure, causing swelling of the biofilm due to 
water intake from the agar gel [36]. Cells are pushed 
from the bulk outwards, speeding up horizontal growth 
and increasing nutrient intake from the agar surface. 

The origin of the wrinkles seems to be controlled by 
the mechanical properties of the biofilm, which are in 
part governed by the production of extracellular matrix. 
In early stages of the evolution of the central biofilm 
core, localized cell death takes place at high density re¬ 
gions due to biochemical stress. Cell death is usually 
followed by cell shrinkage and resorption. Lateral me¬ 
chanical stresses built up during growth are relieved at 
those sites through vertical biofilm buckling, leading to 
wrinkle formation DU- Wrinkles are filled with water 
extracted from agar, as in Figure [3} forming a network 
of channels that enhances nutrient and waste transport 
driven by surface evaporation [37]. The emergence of the 
outer radial wrinkles shown in Fig. [lja) or successive 
wrinkle branching in Fig. [2ja) is less understood. 

Recent developments in the study of the wrinkling of 
gels and thin films on viscoelastic substrata might have 
a counterpart in biofilm growth. A difference with cel¬ 
lular films though, is that there is no feedback between 
mechanical stresses, rates of swelling and growth and cel¬ 
lular differentiation. Swelling of the outer biofilm corona 
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might produce radial wrinkles provided the center of the 
biofilm was stiffer and swelled less than the corona. This 
“corona instability” is commonly observed in swelling 
gels when they are radially graded j39[. Such gels do 
not grow but swell by absorbing water while their elastic 
modulus diminishes. Contact of a film with a viscoelas¬ 
tic substratum may also induce wrinkles, and select the 
dynamics that form them mm- Wrinkle patterns that 
switch smoothly from an intricate core to branching ra¬ 
dial wrinkles have been observed in thin films on a vis¬ 
coelastic substratum due to solvent diffusion in the film, 
which creates residual stresses and stiffness gradients m- 

Based on the previous observations we propose a hy¬ 
brid biofilm description, combining a stochastic treat¬ 
ment of cellular processes coupled to the continuum con¬ 
centration fields and a continuum description of defor¬ 
mation mechanisms informed by cellular activities and 
the interaction with the substratum. Biofilm growth and 
physiology rely on the diffusion and transport of nutri¬ 
ents, waste, and signaling molecules, facilitated by wa¬ 
ter. The elastic deformation of the biofilm is described 
by Foppl-Von Karman equations with residual stresses 
generated by growth and swelling [431145] , formulated for 
a thin biofilm expanding on a substratum 40, 4]| con¬ 
taining water and nutrients. We introduce a strategy to 
estimate residual stresses, as well as spatial variations 
in the elastic parameters of the biofilm, averaging and 
smoothing stochastic information on cell type distribu¬ 
tion, cell division and water absorption. We are able 
to reproduce qualitatively patterns and trends observed 
in experiments, gaining insight on the way biofilms are 
shaped. We argue that stiffness gradients due to growth 
and swelling may cause wrinkle branching and wrinkled 
coronae. We identify governing parameters and study 
their influence on the wrinkled structure. The substra¬ 
tum properties should be carefully measured since they 
largely affect the time scales for biofilm deformation. 
Variations in the values of its viscosity and Poisson ratio 
may change the wrinkling time scales from minutes to 
days, as argued in Section |IV A| Experimental observa¬ 
tions show time variations too, compare Figures 1(a) in 
Refs. [14] and [37j after four and one days, respectively. 
The influence of the agar substratum on the dynamics 
of biofilms and bacterial communities has already been 
noticed in different contexts. Ref. [46] shows the rele¬ 
vance of the agarose concentration in the transition from 
planar to three dimensional behavior during the invasion 
of a three dimensional agar matrix by E. coli. Ref. m 
observes that the water content of the agar substratum 
affects the motility of Pseudomonas strains. 

Section [IT] describes the basic structure of the hybrid 
model. Cellular processes are considered in Section m 
Stochastic division, death, differentiation and EPS pro¬ 
duction mechanisms are introduced. Section HVl discusses 
film deformation according to Foppl-Von Karman equa¬ 
tions and presents water absorption strategies, together 
with procedures to estimate residual stresses and spatial 
variations in the elastic constants. The coupling between 


stochastic and continuous descriptions is detailed in Sec¬ 
tion [V] Simulations of biofilm evolution activating single 
or a few mechanisms are presented throughout the pa¬ 
per. Section [Vl| explores the coupling between processes 
through a few selected tests. Finally, Section |VII| sum¬ 
marizes our conclusions. 


II. HYBRID MODEL 

To investigate the influence of processes at bacterium 
level, cells (or small clusters of similar cells) are repre¬ 
sented by discrete ‘agents’. Different modeling frame¬ 
works describe those ‘agents’ and their evolution in di¬ 
verse ways. Cellular automata distribute them over the 
sites of a regular lattice [48;. The ‘agents’ have inter¬ 
nal state variables, that change according to a set of 
rules, stochastic or deterministic. Standard cellular au¬ 
tomata for biofilms in flows and in porous media allocate 
biomass (cells and EPS), fluids and solids over the lattice 
sites [32] 33 . The dynamics of the biomass in each tile 
is governed by rate equations coupled to concentrations 
and fluid equations. Whenever a threshold is surpassed, 
a neighboring tile is occupied, pushing biomass further 
in the direction of minimal mechanical resistance. Be¬ 
low a threshold, dead occurs with a certain probability. 
Biomass fragments can detach due to erosion by the flow. 
Refs. 33 , 49! use this approach to predict local diffusion 
coefficients, permeability coefficients and Young modulus 
variations taking into account the microstructure. 

Cellular Potts approaches (in particular, the Glazier- 
Graner-Hogeweg variant) have been mostly applied to 
describe aggregation and migration of cells and tissue 
growth, see Ref. m for a review of applications. These 
models use energy based Monte-Carlo updating and are 
able to account for cell volume and shape. This is impor¬ 
tant when describing interactions strongly dependent on 
cell geometry [48]. Effective (not physical) energies typ¬ 
ically include terms representing contact between neigh¬ 
boring cells, a volume constraint on cells, plus chemo- 
taxis terms m ■ Liquids and other substances are rep¬ 
resented as additional cells. According to Ref. [48], the 
DAH (differential adhesion hypothesis) guarantees that 
final shapes correspond to the minima of effective ener¬ 
gies. This hypothesis assumes that an aggregate of cells 
behaves as a mixture of immiscible fluids [48]. 

Individual-based models developed for biofilms in 
aquatic media consider cells as particles (usually spheres) 
distributed in a three dimensional space, whose volume 
and mass grow due to nutrient consumption, see Refs. 
[35] [51] for a survey of applications. Cell division and 
death occur when specified size thresholds are reached. 
Daughter cells are distributed randomly and dead cells 
are removed. The EPS matrix may be represented as 
a fluid continuum, particulates and/or capsular EPS. 
Biomass is eroded according to a law depending on shear. 
Agent growth, division, death and shrinking cause agent 
movement, via shoving algorithms and a pressure field. 
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FIG. 3. (Color online) Schematic representation of a biofilm, 
(a) Biofilm without water network, (b) Biofilm with water 
channels. For computational purposes, space is partitioned 
in a grid of tiles of size a, a being the average length of a 
bacterium. Each tile may be filled with air, water, agar, or 
biofilm biomass (either a bacterium or ECM). This grid is 
used to discretize the continuum concentration and displace¬ 
ment fields, and to track the status of individual bacteria, 
which constitute a fraction of the total biofilm biomass. 


Biomass is advected following a Darcy’s law involving the 
pressure and the dynamic viscosity of the biofilm, which 
is therefore assimilated to a fluid [35|. To investigate the 
effect of the shape of bacteria on the alignment of mi¬ 
crobial communities expanding on a surface, rod-spring 
and sphere-spring representations have been recently im¬ 
plemented m- Each bacterial rod or sphere is charac¬ 
terized by its mass, velocity and position. The EPS ma¬ 
trix is represented by springs that join the rods together, 
and to the substratum. Cells divide when they reach a 
threshold mass due to nutrient consumption. New cells 
and springs must be created in the network. Then, the 
system relaxes to a new equilibrium following Newton’s 
laws. A similar strategy is applied in Ref. [46] to describe 
the influence of the shape of bacteria in the invasion of a 
three dimensional agar matrix by E. coli. 

We are interested here in understanding the elastic de¬ 
formations caused by cell division, death, EPS produc¬ 
tion, swelling, and the interaction with the substratum. 
A framework facilitating the coupling of the relevant mi¬ 
croscopic and macroscopic phenomena is advisable. We 
resort to a cellular automata based hybrid description. 
Hybrid models combine stochastic representations of cel¬ 
lular processes with continuous descriptions of other rele¬ 
vant fields, in our case, concentrations and deformations. 
This allows to take into account the inherent randomness 
of individual bacterial behaviors observed in the exper¬ 
iments uses]. It also permits to consider local inter¬ 
actions around each cell to describe differentiation pro¬ 
cesses. 

In this model, cells are regarded as creatures living in 
a grid, that reproduce, differentiate or die with certain 
probability, informed by the status of the concentration 
fields nearby. In a biofilm, cells are immersed in extracel¬ 
lular matrix. The dominating mass transport mechanism 


in the system is the diffusion produced by the presence 
of chemical concentration gradients of the chemical com¬ 
pounds (oxygen, nutrients, waste products and autoin¬ 
ducers) between the biofilm and the other phases (agar, 
water and air). Cellular processes inform inner mechan¬ 
ical stresses that shape the biofilm. 

We consider a B. subtilis biofilm patch (subdomain 
in Figure [3](a) ) attached to an agar-gel surface (subdo¬ 
main fli) and exposed to the atmosphere (subdomain 
f^)- Eventually, the biofilm may contain channels car¬ 
rying water (subdomain f^) dragged from the agar sub¬ 
stratum, as in Figure ^b). To perform their metabolic 
activities, the biofilm obtains nutrients and other sub¬ 
stances from the gel-water substratum, and from the at¬ 
mosphere. All the subdomains are divided in a grid of 
cubic tiles. The spatial step a chosen for this description 
was the size of one bacterium (around 3-5 fi m) to be able 
to consider individual interactions among bacteria. Each 
tile can be filled with either agar, air, fluid or biofilm 
biomass (either cells or ECM). Partition in tiles allows 
to track the status of individual bacteria, which repre¬ 
sent a fraction of the total biomass, and discretize the 
equations governing the evolution of continuum fields. 



Compute concentration of nutrients/oxygen 
Compute concentration of waste 
Compute concentration of ComX 
Compute concentration of Surfactin 
Compute concentration of EPS 

Detect dead cells 

Detect inert cells Allocate 

Detect dividing cells —► newborn 

Allocate EPS tiles cells 

Update surfactin producers 
Update EPS producers 

Allocate water tiles 
Evaluate liquid flow 

Compute growth tensor 
Compute residual stresses and strains 
Compute displacement fields 
Deform the biofilm-substratum system 


FIG. 4. (Color online) Flow chart of the hybrid model. 


The computational evolution of a biofilm in our model 
is schematized in Figure [4] First, we update the cellular 
status, detecting dividing, dead or differentiated cells and 
allocating newborn cells and water. Then we deform the 
biofilm according to its internal stresses and revise the 
concentration fields and the cellular status. The next 
sections describe how we take into account cellular be¬ 
havior and mechanical processes. 
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III. DESCRIPTION OF CELLULAR 
PROCESSES 

It is commonly accepted that biofilms are able to gen¬ 
erate molecular signals as cell-cell communication mech¬ 
anisms to activate gene expressions that unleash differ¬ 
ent survival strategies adapted to the environmental con¬ 
ditions mm- The basic approach relies on the idea 
that the same organism that is sensible to certain chem¬ 
ical inducers is also responsible of producing those com¬ 
pounds, creating a self-feeding loop which regulates both 
the genetic expression and the generation of autoinducer. 
Quorum sensing mechanisms provide a way for bacteria 
to count their population prior to unleashing different 
behaviors such as differentiation, virulence or spreading 

[EM- 

Differentiation is a mechanism through which some 
bacteria inside a population display physiological alter¬ 
ations, developing different phenotypes without a change 
in their genotype. These modifications lead to the de¬ 
velopment of specialized cells with certain characteristics 
which perform specific tasks. 


A 



SURFACTIN EPS CELLS 


Nutrient /Oxygen 
shortage 



NORMAL INERT 

CELLS CELLS 


FIG. 5. (Color online) Directional signaling in B. subtilis. 
When a threshold concentration of ComX is reached, a per¬ 
centage of the population of cells expresses genes that induce 
a change in their phenotype, becoming surfactin producers. 
Surfactin may signal other cells to produce extracellular ma¬ 
trix. Starving cells become inert spores. 

Differentiation processes have been observed to depend 
on the local conditions around each cell, in particular, the 
state of the neighbours and the local chemical composi¬ 
tion found around them mu [57]. B. subtilis provides a 
well known example, where initial biofilm cells may differ¬ 
entiate in four well defined cell types. Two autoinducers, 
ComX and surfactin, relate the four types of differenti¬ 
ated cells [14], see Figure [5j The different cell types are: 
• Normal cells. They perform all metabolic activities 


and generate ComX. They are sensible to all autoinduc¬ 
ers and may be deactivated due to a lack of nutrients or 
oxygen. They may turn into surfactin generators or cells 
expressing matrix genes when threshold concentrations 
of ComX or surfactin are reached. 

• Surfactin producers. They generate surfactin and do 
not reproduce. Proliferation of these cells may be inhib¬ 
ited by the presence of EPS producers in their neighbor¬ 
hood. 

• EPS matrix producers. They secrete EPS matrix, and 
are able to reproduce at lower rate than normal cells 
because their metabolism also spends resources on gen¬ 
erating EPS. 

• Inert cells, or spores. These cells do not perform any 
metabolic activities, neither reproduction, nor differen¬ 
tiation. EPS producers become spores when there is a 
severe lack of nutrients in the local environment of the 
cell. Spores may activate again when the environmental 
conditions improve. 

As time progresses, a fraction of cells expresses surfactin 
genes. Later on, an increasing number of cells secretes 
EPS matrix. Cells may also die as a result of biochem¬ 
ical stresses m- This has been shown to occur at the 
interface agar-biofilm as a result of biochemical stresses 
at localized high density regions. Eventually, sporulation 
starts and few normal cells remain. Spores are usually 
located in upper regions of the biofilm, whereas normal 
cells are restricted to the bottom edges m- 

We describe below the rules for cell differentiation, 
death and reproduction. 


A. Cellular division, deactivation and death 

When there is an excess of oxygen, the concentration 
of nutrients c n becomes the limiting concentration that 
restricts biofilm growth inside the biofilm D 3 : 

Cn,t D n i)N.C n = k n ■ -zrz—. ( 1 ) 

C n + K n 

Nutrient uptake is simply described by a Monod law. 
Inside the agar substratum Oi: 

Gi,t -^n,a^C n = 0- (2) 

In the water phase D 4 (if present): 

C n ,t - D n jAc n + V • Vc n = 0. (3) 

Here, D n ^, D n a and D n ^ are diffusion constants in the 
biofilm, agar and water, respectively, K n is the Monod 
half-saturation coefficient, and k n {x) is the uptake rate of 
the nutrient at the particular location. It will be set equal 
to k n in tiles occupied by alive cells and equal to zero 
elsewhere, v is the transport velocity in the water chan¬ 
nels. On the interfaces agar-biofilm and water-biofilm 
(if present) we impose transmission boundary conditions 
(continuity of concentrations and fluxes). On the inter¬ 
face biofilm-air, the interface agar-air and the edges of 
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the agar substratum, we have zero flux conditions. Ini¬ 
tially, the nutrient concentration is usually taken to be 
constant c n = C in agar. 

Tiles occupied by alive cells C are assumed to divide 
with probability ESQ: 


Pd(C) 


c n (C) 

c n{C) + K n 


(4) 


c n being the limiting concentration. Initially, we real¬ 
locate newborn cells inside biofilm by pushing existing 
cells in the direction of minimum mechanical resistance, 
that is, the shortest distance to the biofilm-air interface 

mm- 

Figure [ 6 ] illustrates the expansion of a biofilm seed. 
The initial round cluster grows in size in a radial way. 
Brick-like colonies are formed. Their border may advance 
in a wavy way depending on the velocity, as usual when a 
free boundary expands outwards [58] . This spread mech¬ 
anism may change as the biofilm develops and water is 
absorbed, as explained later. 


(a) (b) 



.0 


FIG. 6. (Color online) Expansion of a circular biofilm seed 
driven only by nutrient consumption. The initial seed is 
formed by two layers of cells with a diameter of 40 cells. 
Only division and spread mechanisms are involved. Snapshots 
taken every 100 steps. The initial nutrient concentration in 

h 2 

the substratum is C — 3 K n and the ratio D nC ^ — 8. 


The constants K n affecting the probability laws are 
numeric constants that should be calibrated with experi¬ 
mental data to have a real description in the simulations 
and do not have to match with saturation constants used 
in the source terms of the react ion-diffusion equations for 
the concentrations, which are related with a chemical sat¬ 
uration limit. To simplify, they have been chosen to agree 
with the values for saturation constants. 


Cells may suffer a shortage of nutrients, triggering a 
sporulation mechanism to preserve bacteria in a “deacti¬ 
vated” state until the environmental conditions improve. 
To reflect this, cells with a local limiting concentration 
below a certain threshold will become spores in our simu¬ 
lations. We deactivate cells C with probability 1 — Pd(C) 
whenever Pd(C) < 7 d for a certain threshold 7 ^ > 0 . 
These cells may activate again when the conditions im¬ 
prove. 

As mentioned before, cells may also die due to bio¬ 
chemical stress in high density regions where waste prod¬ 
ucts and toxins accumulate but nutrients deplete. Taking 
for instance the concentration of waste c w at a location 
as an indicator of death, a cell C is scheduled to die with 
probability: 


P W (C) = 


AC) 


Cw(C) + K w ’ 


(5) 


whenever P W (C) > 7 w for a certain threshold 7 ^ > 0 . 
The concentration of waste is governed by 


u,t Dw,b^Cw k w ^ ^ 


(6) 


inside the biofilm ^ 3 . Si >w (x) equals 1 in regions occupied 
by alive cells C{ and vanishes otherwise. On agar-biofilm 
and air-biofilm interfaces we impose no flux boundary 
conditions. If a fluid phase ^4 is present, we impose 
transmission boundary conditions at the interface and 
include a diffusion-convection equation similar to Eq. © 
for waste transport in the liquid. D w D w j are the 
diffusion constants, and k w the production rate of waste 
and toxins at the particular location. 


B. Surfactin and EPS production 

Normal cells may turn into surfactin generators when 
a threshold concentration of ComX is reached. ComX 
is a molecule related with quorum sensing mechanisms 
in B. subtilis [14] . It is produced by all the cells in the 
biofilm, except those which are not metabolically active. 
The concentration c cx of ComX inside the biofilm ^3 is 
governed by: 

£cx,t b Ac cx &CX ( 1 j ^ ^ CX5 (7) 

V C cx T -t\-cx / 

l 

with no flux boundary condition on agar-biofilm and air- 
biofilm interfaces. & cx is the production rate, and D C x,b 
the diffusivity. Si, cx (x) equals 1 if cell produces ComX 
and vanishes otherwise. k cx is multiplied by an inhibition 
factor for large concentrations. A guess for K cx may be 
the half-saturation constant of ComX. If a liquid phase 
f ^4 is present we impose transmission boundary condi¬ 
tions at the interface and include a convection-diffusion 
equation similar to Eq. © for ComX transport in the 
liquid. Normal cells C are assumed to become surfactin 






producers with probability: 


Ps(C) = 


Ccx(C) 

Ccx(C) + i^cx ’ 


( 8 ) 


whenever P S (C) > 7 S . The threshold y s reflects that dif¬ 
ferentiation should start when a minimum background 
concentration of ComX is reached. 

When EPS producers are present, they act as in¬ 
hibitors. We may set P s = P in hib Cc °+k cx with p inhib 
the ratio of neighbors producing EP§ to total number 
of neighbors. If a cell is surrounded by EPS producers, 
there is no any chance of differentiation into a surfactin 
generator. A similar effect is obtained setting a much 
lower threshold for differentiation into a surfactin pro¬ 
ducer compared to EPS producers, as we usually do in 
our simulations. 

Surfactin acts as an autoinducer in normal bacteria, 
changing their phenotype to become EPS producers [36] . 
The concentration c s of surfactin inside the biofilm is 
governed by: 


G ,t 


Ds,b^Cs — k 



(9) 


with no flux boundary condition on agar-biofilm and air- 
biofilm interfaces. k s is the production rate, D s ^ the dif- 
fusivity and K s the half-saturation of surfactin. ^ ?s (x) 
equals 1 if cell Ci produces surfactin and vanishes oth¬ 
erwise. When a liquid phase Q 4 is present transmis¬ 
sion boundary conditions hold at the interface and a 
convection-diffusion equation similar to Eq. ©> governs 
surfactin transport in the liquid. We assume that ac¬ 
tive cells C not secreting surfactin become EPS producers 
with probability: 


Pe(C) = 



.(C) 


Cn(^) T Kn 


( c s (C) \ 
\c s (C) + K s ) 5 


( 10 ) 


whenever c (q+k > 7ei> which means that a minimum 
concentration of surfactin is required. The probability 
increases with the availability of surfactin and the scarce¬ 
ness of nutrients. This is consistent with measurements 
reported in Ref. [59]. If we want to enforce a level of 
nutrient depletion to allow differentiation into EPS pro¬ 
ducers we may further impose c (c)+k < 7e 2 • This may 
forbid differentiation below a certain height. 

Figure [7] illustrates the distribution of cells when we ac¬ 
tivate the differentiation processes, in combination with 
standard growth and spread. As the size of the biofilm 
increases, the concentration of ComX raises until the 
threshold is reached. At this point bacteria start to 
differentiate and secrete surfactin. This triggers differ¬ 
entiation of cells into EPS producers when the surfactin 
threshold y ei is reached. As the number of EPS cells 
increases, the rate of differentiation into surfactin gen¬ 
erators stabilizes and may even decay due to an inhibi¬ 
tion mechanism induced by the EPS cells M- Relative 
proportions for each type of cell were adjusted by fitting 


10 cells 
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FIG. 7. (Color online) Slices of a biofilm showing the cell type 
distribution during early stages of the differentiation process. 
All the slices correspond to the same biofilm and are taken 
at the same time. The evolution started from an initial seed 
with a diameter of 60 cells and variable thickness, represented 
by three small peaks. Dead cells appear first where the initial 
peaks where located. In later stages, inert cells proliferate at 
the top whereas normal cells are restricted to the bottom and 
the edges. Each colored box is one cell. Parameter values: 
C = K n , = 0.01, = 0.001, = 0.01, 


D w K vl 


k s a z 


D K — 0.8, 7 ™ — 0.005, 7 d — 0.00001, - 0.001, 7 ei — 

0.60601,7e 2 = 0.2. 


parameters, but there should be an experimental calibra¬ 
tion to improve these values. Normal cells are found in 
the edges and near the interface biofilm-agar. EPS pro¬ 
ducers appear once surfactin is generated and a certain 
thickness is reached so that the nutrient concentration 
depletes. Dead cells concentrate on the bottom surface 
at inhomogeneity peaks when the initial seed has variable 
height, due to higher waste concentration. As the biofilm 
increases its thickness, scarcity of nutrients is felt in the 
upper part due to the high resistance to mass transfer 
across the increasing thickness. This effect induces the 
appearance of spores in the upper part, slowing down the 
growth of the biofilm. The geometric distribution of each 
type of cell seems to follow qualitatively trends reported 

in Refs. QUEUES]- 

Implementing the differentiation strategy described 
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above requires measuring a number of kinetic constants 
for different concentrations. Since these values are cur¬ 
rently unknown, we have scaled the parameters to obtain 
the expected behaviors in small clusters of cells, lower¬ 
ing the computational cost. For Gram negative Vibrio 
harveyi , quorum sensing processes triggering biolumines¬ 
cence have been modeled in terms of feedback loops (56] 
and applying information theory to the quorum sensing 
circuit m, which may reduce the number of parameters 
involved. Signaling parameters have been deduced from 
in vivo data for V. harveyi [61] . 

For simplicity, the biofilm in Figure [7] does not include 
the ECM phase, only cells. We might include it following 
the procedure described below. Once EPS production 
starts, a new concentration c e diffuses in the environ¬ 
ment: 


Ce,t De,b^C e — k e 



C e 

Ce + K e 



(ii) 


k e is the production rate, D e ^ the diffusivity and K e the 
half-saturation of EPS. ^ ?e (x) equals 1 if cell C; is an EPS 
producer and vanishes otherwise. Initially, these sub¬ 
stances are just compounds diffusing in the environment. 
As the concentration increases they become a polymeric 
envelop for the cells. This description can be improved 
taking into account the presence of molecules of different 
size. Monomers diffuse, but longer molecules stay where 
they are produced. The concentration of momoners is 
governed by a diffusion equation coupled to the concen¬ 
tration of larger molecules, whereas the latter is governed 
by rate equations, that substitute Eq. ©• 

The volume fraction of biofilm occupied by extracel¬ 
lular substances varies largely depending on the type 
of biofilm. Images in Ref. [13 j show biofilms in flows, 
formed by scattered cells embedded in a dominant EPS 
phase. On the contrary, imaging biofilms on air-semisolid 
interfaces at hight resolution shows a structure of densely 
packed cells glued by EPS, see references m, m and 
Ej. Space is mostly occupied by cells. This suggests 
skipping the creation of an EPS phase and thinking of it 
as a virtual glue that keeps cells together and modifies 
the mechanical properties of the biofilm, in particular, its 
elastic constants and the osmotic pressure. If we choose 
to create an EPS phase, we may exploit the concentra¬ 
tion c e of EPS to generate new biofilm tiles filled with it. 
At any tile, we set 


Pm(C) 


Ce(C) 

C e (C) + K e 


( 12 ) 


If this value is larger that a threshold y m , we create an 
EPS tile at that location and shift biomass around in the 
direction of least mechanical resistance. should be 
calibrated so that the fractions of biomass corresponding 
to cells and EPS fit the experimental observations. The 
true cells are a fraction of the total amount of tiles occu¬ 
pied by biomass. They divide and differentiate according 
to the same probabilistic rules we have already described. 


Individual-based models for biofilms on agar surfaces 
represent the EPS as springs joining the cells [51]. For 
biofilms in flows, expressing larger fractions of EPS, each 
‘individual’ contains cellular mass and capsular EPS, 
that vary according to user defined rate equations. Ad¬ 
ditional EPS particulates may be created as EPS pro¬ 
duction increases following user defined rules [35] . Refs. 
[221162] study EPS dynamics in a quorum sensing frame¬ 
work, defining competing strains: EPS producers, EPS 
non producers, strains that downregulate EPS produc¬ 
tion for large cell densities. This is not the reported quo¬ 
rum sensing circuit for B. subtilis , see Ref. [14] . Cellu¬ 
lar automata models for biofilms in flows distribute cel¬ 
lular mass and EPS in each lattice site following phe¬ 
nomenological rate equations (32, 33]. Spread to neigh¬ 
boring sites occurs as a threshold mass is surpassed. In 
both frameworks, cells are scheduled to die when their 
mass or volume falls below a threshold. No memory of 
their former presence, relevant for elastic deformations, 
is kept. Neither the relevance of waste accumulation for 
cell death, pointed out in Ref. m, nor the currently 
known quorum sensing circuits for EPS production are 
incorporated in those models. Dynamic energy budget 
descriptions of cell dynamics [63] might provide a sophis- 
phicated framework to describe the evolution of individ¬ 
ual cells taking into account their mass, volume, mainte¬ 
nance, substance secretion and the presence of nutrients, 
oxygen and toxicants. However, extensions to biofilms 
are yet to be devised. 


C. Nondimensionalization of the concentration 
equations 

Nondimensionalization of the equations for the con¬ 
centrations listed above reveals a key governing param- 

7 2 

eter for each concentration: F = where k is the 

uptake or production rate, a the bacterium size, K the 
half-saturation constant and D 5 the diffusivity inside the 
biofilm. Setting c = and x = the equations for 
concentrations become: 

(13) 

inside the biofilm, where / represents the production or 
consumption terms, and 

- AxC + v • VxC = 0, (14) 

in the water channels, where v = 7=^. We have sup- 
pressed the time dependence because diffusion of chemi¬ 
cal species is faster than the cellular processes considered. 
For each biofilm configuration, we solve these nondimen- 
sionalized stationary equations with the corresponding 
boundary conditions by a relaxation technique. Then 
we compute all the probabilities for division, death and 
differentiation, detect the different types of cells and allo¬ 
cate newborn cells. This generates a new biofilm config¬ 
uration for which the procedure is repeated. A complete 
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reproductive cycle of a bacterium was set as the basic 
time step for the cellular processes in the model (which 
is typically in the order of 26-45 minutes depending on 
nutrient and temperature). This choice implies that our 
system is discrete in space and time [SD- 


iv. DESCRIPTION OF MECHANICAL 
PROCESSES 

Differentiated cells produce chemicals, which may in¬ 
duce further differentiation of other cells or trigger me¬ 
chanical processes that alter the biofilm structure, im¬ 
proving its chances of spread and survival. EPS produc¬ 
tion increases the stiffness of the biofilm m, which may 
be identified with a material undergoing macroscopic de¬ 
formations caused by growth [16] over a substratum. The 
spatial distribution of dead and inert cells, together with 
water absorption, may induce additional stiffness gradi¬ 
ents. This motivates allowing our biofilm to move along 
the tiles according to microscopically informed contin¬ 
uum descriptions of the underlying mechanical processes. 


A. Out-of-plane deformations and wrinkle 
formation 


Due to the interaction of components in the extracellu¬ 
lar matrix, dissipative phenomena take place inside soft 
tissues during deformation processes and the mechani¬ 
cal response of tissues if often viscoelastic. Nevertheless, 
growth induced deformations on time scales much larger 
than those of relaxation may be regarded as purely elas¬ 
tic [45]. Since the height of real biofilms is about 10-100 
times smaller than their radius and they are initially thin 
and flat, we may use Foppl-von Karman [43] equations 
for thin elastic plates to study their elastic deformation 
reducing dimensionality [ljy. Residual stresses due to 
growth are incorporated in the description in terms of 
the growth tensor, see Refs. [44j[45]. The shape of the 
body is then determined by both growth and elastic de¬ 
formation. The equilibrium deformation of a growing tis¬ 
sue is described by a variant of the Foppl-Von Karman 
equations [45] : 


D (AH - A C M ) - = p, 

dxp ~ ’ 


(15) 

(16) 


where a,/? stand for x,y and summation over repeated 
indexes is intended. The first equation describes out-of- 
plane bending £(x,y) and the second one in-plane stretch¬ 
ing for the displacements u = (u x (x,y),u y (x,y)). (x,y) 
vary along the 2D projection of the 3D biofilm structure. 
The bending stiffness is Dm ? anc ^ ^ initial 

plate height. P is the external pressure at the border of 
the sample. Cm = d( ^ 9zx ^ 9xz ^ + d ^z y +^g yz ) - g a res iq ua i 


growth term. The components of the tensor g are the 
growth rates, that is, the rate of volume supply per unit 
volume. Defining consistently a growth tensor g is a non 
trivial issue. A possibility is to set g = Vw m, w being 
the vector field of displacement of growing biomass in the 
biofilm w, which is determined by the cellular division 
and spread mechanism and water absorption processes. 

Stresses a and strains 6 are defined in terms of in-plane 
displacements u = (u x ,u y ) [39[ [45] by: 


z cx.,(3 


OXX - 


1 (du a du 0 d£ d£ 


+ ^ + 


+ e a,/3> 


2 \dxp dx a dx a dxp 
E ^ _ _ - 

l _ v £xx E V£ yy)i (Jx y ~ \j rV e - 


E 


xy> 


E 


J yy 


1 - z /2 ^ yy 


(£yy T VE X x)' 


(17) 

(18) 


The Poisson ratio of different tissues has been measured 
to be about 0.49999 [65] . We will set the Poisson ratio of 
the biofilm v = | (incompressibility) f45]. The average 
elastic modulus E of a wild type B. subtilis biofilm has 
been measured to be about 25 kPa P3- When the resid¬ 
ual strains a are due to growth, they are expressed in 
terms of the growth tensor as: 

£ a,(3 = ~2 (y a/3 + y /3a + 9zcx9zg) • (19) 

A growing biofilm is in a state of compression due to cell 
division and, eventually, water absorption. Alternatively, 
this may be represented by a residual strain 


£ ° a ,p( x ^) = -£o(x-,t)S a ,p, £o > 0. (20) 


If we assume that cells do not grow at expense of their 
neighbors, 5 a ^p is a diagonal unit tensor in polar coordi¬ 
nates in a circular film. 

For a biofilm growing on a surface no external loads act 
on the edges, therefore, P = 0. However, the interaction 
with the surface affects the deformation process and has 
to be included in the description. A correction to Eqs. 
(15)-(16) applicable to thin elastic films growing on a 
viscoelastic substratum is proposed in Ref. m : 


dt 


1 — 2 v v h v 

2(1 z/i>) 9v 


D(- A 2 £ + A C M ) 


+ h 


d 

dxp 



du 

dt 


h v h 

r\v 


V • <t(u) - —u, 



( 21 ) 

( 22 ) 


where h v is the thickness of the viscoelastic substratum 
and fjb v , v v , g v its rubbery modulus, Poisson ratio, and 
viscosity, respectively. 

The equilibrium equations describe possible equilib¬ 
rium configurations, and yield information on their sta¬ 
bility or changes of stability. To describe the dynam¬ 
ics of wrinkle formation and the finally selected pattern 
we must solve the time dependent equations. Equations 
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( 21 )-( 22 ) were derived for the interface between the thin 


film and the substratum. This interface moves vertically 
following the displacement field £. This shifts all the tiles 
in the three dimensional biofilm-substratum system ver¬ 
tically. 

Choosing new dimensionless variables: 

x = z- u = z- ( =h’ ° = E’ t = f■ 



and setting L = 7 / 1 , the dimensionless equations read: 


d[ 

dt 


12(1 - i/ 2 )7 2 


d 

dx 0 



du 

dt 


(-A?f + A xCm) 


- 

Vv 


T Vx • ct(u) - T — u, 
Vv 


with 


T = 2(1 - v v ) rjyh 12(1 - z/ 2 ) 7 4 
1 — 2v v h v E 


(23) 

(24) 


(25) 


We set L equal to the maximum averaged radius of the 
biofilm in study. Common experimental values for the 
parameters defining the spatial scales are h ~ 100 /im, 
h v ~ 100 h. The size of the biofilm may vary with the 
nutrients and the agar gel nature. Figure 1(a) in Ref. 
m shows a wrinkled biofilm with L ~ 50 h grown in 
four days on a 1.5% agar surface. Figure 1 (a) in Ref. 
m reproduces wrinkled biofilms with L ~ 50 h within 
one day. 

The linear terms T ^-u and T are just damping 


contributions from the substratum slowing down the evo¬ 
lution. We will neglect them in our preliminary stud¬ 
ies. Wrinkling behavior is observed when the coefficient 
12(1 — v 2 )^ 2 is large enough, so that the contribution 
of nonlinear effects and residual stresses in Eq. (23) be¬ 
comes relevant and drives the plate out of the flat equi¬ 
librium state, see Figure | 8 | 

The factor r in Eq. (24) separates the time scales for 
the evolution of u and £. The closer to 0.5 the Poisson 
ratio v v is, the larger this separation is. The in-plane dis¬ 
placements u may reach a quasistatic situation in which 
their evolution is driven by slower changes in the resid¬ 
ual stresses created by growth, swelling, and off-plane 
displacements. 

The constitutive parameters for an agar gel may vary 
depending on the agar source, its content of agar and 
other products, in particular, water. When an hydro¬ 
gel is fully swollen, its mechanical behavior is similar to 
those of rubber-like materials [ 66 ] . which have a Poisson 
ratio of about 0.5. For some agar gelatines the Young 
modulus is measured to be 0.4999, with an elastic mod¬ 
ulus about 27 kPa [65]. In other agar gels these numbers 


FIG. 8 . (Color online) Connected network of wrinkles switch¬ 
ing from an intricate core to radial branches, formed in a ho¬ 
mogeneous circular film when a circular front of residual com¬ 
pression stresses of magnitude £0 = 0.1 expands at constant 
velocity. Parameter values: v — 0.5, E — 25 kPa , v v — 0.45, 
/i v = 0 and 7 = 16. dx = 0.1 h , where h is the film thickness 
before wrinkling. 


go down to v v = 0.32 with an elastic modulus E v = 52 
kPa [BY . This modulus increases with the agar concen¬ 
tration from ~ 30 kPa in 0.5% samples up to 700 kPa in 
5% samples. The viscosity r] v of agar gels is reported in 
the range of 100 cp = 0.1 Pa-s for 1.5% samples, see Ref. 
[ 68 ]. As the agarose concentration increases to 8 %, the 
viscosity raises to 1 — 1.5 Pa-s [69]. These values increase 
as the temperature diminishes. According to Ref. 69l . 
no measurements can be obtained below 36° C. Standard 
experiments are carried out at room temperature, usu¬ 
ally below that threshold. Effective values for viscosities 
in a MPa-s range are adjusted in Ref. m • 

A large variability in the time scales may arise due 
to uncertainty in the experimental values of v v and r] v . 
Setting for instance v v — 0.4999, r\ v =0.1 Pa-s, we find 
T ~ 31 hours. This value varies enormously depending 
on that changes with the water content. Switching to 
v v = 0.45 and rj v = 1 Pa-s, we find T ~ 2480 s ^ 40 m. 
The time scale scale for u is T' = ^ ranges from 10 -3 s 
in the first case to 10 -2 s in the second one. If r\ v enters 
the range of kPa-s or MPa-s, this value increases by a 
factor 10 3 — 10 6 . 


B. Water absorption and water channels 

As mentioned earlier, variations in the osmotic pres¬ 
sure extract water from the agar substratum. In Ref. 
[3B], a macroscopic model for water migration from the 
agar substratum into the biofilm is proposed in terms 
of the volume fractions of water and biofilm at each lo¬ 
cation. The biofilm is formed by cells and extracellular 
matrix, which is often considered a gel with ability to 
swell. As observed in Ref. m, flat regions of the biofilm 
are very resistant to flow. Water cannot be driven inside 
without fracturing the biofilm or delaminating it from 
the substratum. By contrast, water flows beneath the 
wrinkles, forming an intricate network of channels m • 
Wrinkles origin in the core of the biofilm associated with 
dead areas El- Water flows along the channels driven 
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by surface evaporation m- 

These observations suggest a strategy to incorporate 
water into our biofilm description. Water absorption 
by the biofilm increases the volume of the cells and the 
EPS phase. Enlarging the size of the tiles to reflect that 
fact is unpractical. We resort instead to inserting water 
tiles in the biofilm, representing water absorption by the 
biomass. In that way, the biofilm contains a volume frac¬ 
tion of adsorbed water that changes its volume. A low 
cost strategy, easy to implement in our hybrid framework 
is the following. As in Ref. [36] . we propose for the os¬ 
motic pressure a law of the form 7 r = II</>, but replacing 
the volume fraction <p of biomass by our information of 
the biomass available at biofilm columns 

N 

= n—- , ( 26 ) 

-L' max 

where N is the number of tiles occupied by biomass in 
a column and N max is the total number of tiles in the 
column, including water. Water tiles are created with 
probability: 


Pi{C) = 


7 n{C) + E{Cy 


(27) 


where E(C) is the Young modulus at tile C and the height 
of the columns in Eq. (26) is updated as we create wa¬ 
ter tiles. The status of neighboring tiles is shifted in the 
direction of minimal mechanical resistance, except when 
water occupies a dead cell. More water will be absorbed 
in the columns containing a larger fraction of biomass, 
and in softer regions. Scattered inner tiles do not consti¬ 
tute a separate phase, but are considered to be part of 
the swollen biomass. 

We compute E(C) modulating the measured reference 
value of the Young modulus of the biofilm with an av¬ 
erage of weights representing status of the neighboring 
tiles. Weights range from zero for water, a small posi¬ 
tive value for dead cells, a larger value for normal cells, 
up to one for EPS producers. Alive cells attached to the 
substratum contribute much larger weights than the rest. 
E(C ) will be small wherever we have dead cells. As Fig¬ 
ure [9] shows, those regions fill with water in successive 
steps, becoming a water phase that expands following 
the wrinkles. Combining water transport through those 
channels, with the cellular processes described in Section 
m we see that cell death due to waste accumulation is 
reduced and the colony is able to expand in a sustained 
way. In early stages we do not allow newborn cells into 
this water phase due to the pressure gradients that drift 
water from the agar gel into the biofilm. As the water 
content of the agar gel diminishes we might allow cell 
expansion around the channels, as observed in Ref. m 
This would keep water inside the biofilm. 

In Figure [9] we have distributed dead cells above the 
substratum, in the wrinkled regions. The value of II is 
chosen small enough to avoid water proliferation inside 
the biofilm, about E/ 20. These simple probability laws 
produce compression residual stresses that vary with the 
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FIG. 9. (Color online) (a) Biofilm-substratum system de¬ 
formed following out-of-plane displacements £ predicted by 
Eqs. ([ 2 T]) - ([22]) for residual stresses expanding as a radial com¬ 
pression front, (b) Bottom slice of (a) displaying the arrange¬ 
ment of biofilm and substratum, (c) Intemediate slice of (a) 
showing the allocation of dead cells, (d) and (e) Accumula¬ 
tion of absorbed water in slice (c) triggered by the presence 
of weakened dead areas after 1 and 40 steps of the water ab¬ 
sorption process, respectively. 
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where II = E/3. 


biofilm height, as shown in Figure 
Values in the range II ~ E are suggested in Ref. m- 


A time scale for water absorption processes may be 
inferred from the analysis in Ref. 36], where a thin film 
description of biofilm spread due to swelling and growth 
is proposed. For radial biofilms and small departures 
from osmotic balance, a similarity approximation for the 
height h and the radius R of the biofilm yields a time 
scale 1 /q for the spread process, where q is the rate of 
biomass production. 


We cannot directly couple the model in Ref. [36 to a 
cellular automata description of cell division, death and 
differentiation to describe biofilm spread due to swelling. 
The main reason is that it already includes a source rep¬ 
resenting growth. An alternative model of fluid transport 
in the biofilm would be necessary for a better description 
of the water absorption process and biofilm spread in our 
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framework. 

Once water channels have been carved in the biofilm, 
we might couple our cellular automata description to the 
hydrodynamics of the flow as done for biofilms in porous 
media (33]. This assumes that the deformation of the 
biomass and the extracellular fluid flow can be computed 
in a sequential manner, which happens when the biomass 
behaves as an elastic material or the flow induced de¬ 
formation of the biomass matrix is negligible [71]. Us¬ 
ing order-of-magnitude analysis, Ref. m establishes do¬ 
mains of validity of different models for liquid transport 
in a fluid-solid system depending on the volume fraction 
and viscosity of the fluid, the volume fraction, elastic con¬ 
stants and density of the solid, the hydraulic permeability 
of the system, the characteristic times for displacement 
of the solid, and the characteristic macroscopic length of 
the system. 


C. Growth tensor and elastic constants 


The components of the tensor g entering Eqs. (23)-(24) 


are the growth rates, that is, the rate of volume supply 
per unit volume. New tiles may add to our biofilm due to 
cell division, EPS creation or water adsorption. Wher¬ 
ever they are created, they shift another tile in the direc¬ 
tion of minimal mechanical resistance. The growth ten¬ 
sor at each grid location (ix, iy , iz) is computed by keep¬ 
ing track of all the new tiles inserted and the direction 
in which their predecessors where shifted. We define a 
vector w = (wi(ix,iy,iz),W 2 (ix,iy,iz),ws(ix,iy,iz))a. 
wi(ix,iy,iz) is determined by cumulatively adding ±1 
for each tile shifted in the x direction in the positive or 
negative sense, respectively, and w 3 are evaluated in 
a similar way, along the y and 2 directions. The final 
vector is normalized to have norm a. Then, we com¬ 
pute Vw, where the derivatives are approximated by fi¬ 
nite differences that use the known grid values. To es¬ 
timate g(ix : iy) we consider all the contributions from 
Vw (ix,iy,iz) for varying iz. 

The resulting tensor usually presents abrupt oscilla¬ 
tions due to stochasticity, constrained motion along a 
restricted set of directions in a cubic grid and varying 
biofilm thickness. Such oscillations in the grid scale may 
destabilize numerical solutions of Eq. (23). Averaging 


and smoothing these tensors to avoid numerical instabil¬ 
ities, we obtain residual stresses that serve as a basis for 
the numerical tests, see Figure [lO] Averaging rapidly os¬ 
cillating coefficients and sources is a standard practice to 
study macroscopic deformations of materials. 

The reference Young modulus of the biofilm is esti¬ 
mated from experimental measurements PH- The knowl¬ 
edge of its microstructure gained from our simulations al¬ 
lows to introduce spatial modulations. At each location, 
we may regulate the reference value multiplying it by the 
sum of weights describing the status of its neighboring 
tiles, as we did in section |IVB A possible qualitative 
choice for those weights are zero for air, a small fraction 
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FIG. 10. (Color online) Residual stresses s xx f° r the biofilm 
seeds depicted in (a), (b) and (c). (d), (e) and (f) represent 
their values after one step of the stochastic division and spread 

h 2 

process, with C — 2K n and D n ° K — 0.01. Spatial variations 
related to thickness are identifiable in (g), (h), (i), obtained 
averaging a few trials of the stochastic procedure. These fields 
approach their average values (about —0.1). A single trial 
already gives the average value used in rough approximations 
by constants, (j), (k), and (1) show the effect of including dead 
spots in the biofilm seed. Localized depressions are observed, 
(m), (n), (o) reproduce the residual stresses generated by the 
water absorption scheme when II = E/3. The average values 
are about —0.2. s yy resembles e xx in all cases. £ xy is about 
3-2 times smaller. 


of one for dead cells, a larger fraction for normal cells 
and one for cells producing EPS. The biofilm becomes 
softer where it swells easily. In this way, we obtain spa¬ 
tially dependent elastic moduli that allow to investigate 
variations in the macroscopic deformation of the biofilm 
in response to changes in its microstructure. 
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V. COUPLING OF CELLULAR AND 
MECHANICAL PROCESSES 


The interaction between the previous descriptions of 
cellular and mechanical processes is visualized in Figure 
[4j The scheme below details a practical implementation. 


1. Initialization. A matrix is created to indicate the 
status of the tiles in the computational grid of step 
a. At first, we usually set S(ix,iy,iz) = 2 in the 
regions occupied by a biofilm seed, assuming that it 
is formed by undifferentiated cells stuck together. 
As tiles (ix, iy , iz) are filled with differentiated cells, 
EPS or water in the next steps, S will take different 
positive integer values. Sections occupied by air or 
substratum correspond to zero or negative values. 
This grid is used to keep track of cellular processes 
and also to discretize continuum equations. As the 
biofilm enlarges, we might resort to a coarser grid 
for the latter purpose. Initially, the concentration 
of nutrients is set equal to a constant C in the sub¬ 
stratum and zero elsewhere. The concentrations of 
waste, ComX, surfactin and EPS are set equal to 
zero everywhere. 


2 . Cellular processes. 


(a) Concentration update. The concentrations of 
nutrients, waste and ComX are evaluated solving 
Eqs. ([l])-(|3|, (| 6 | and (|7|) with the specified bound¬ 
ary and initial conditions. Whenever surfactin or 
EPS producers are present, Eqs. © and are 
solved too. All the equations are nondimensional- 
ized as indicated in expressions (fl3|)-([l4|). The so¬ 
lutions of the time dependent diffusion problems for 
these chemicals relax to their stationary states in a 
short time, compared to the typical time scales for 
cellular processes. Explicit finite difference schemes 
provide a low cost approximation to the stationary 
concentrations after a number of steps. Those ap¬ 
proximations are stored and used in the next stages 
to evaluate behavioral probabilities. 


(b) Probabilities for cell activities. They are only 
computed at tiles (ix, iy , iz) occupied by cells. Cells 
that are not already dead are killed with probabil¬ 
ity P w whenever the concentration of waste sur¬ 
passes a threshold 7 ^. We evaluate expression 
([5]) and generate a random number r G [0,1]. If 
max(r, 7 ^) < P w , we kill the cell. Alive cells de¬ 
activate with probability 1 — P^, Pd given by for¬ 
mula Q, as long as the nutrient concentration is 
low enough. Active cells become surfactin produc¬ 
ers with probability P s defined in formula (J 8 | if the 
concentration of ComX surpasses a threshold. Ac¬ 
tive cells not releasing surfactin become EPS pro¬ 
ducers with probability P e given by formula (10) 
when a minimum surfactin level is reached and the 
nutrient concentration is depleted. 


(c) Cellular division and spread. Active cells not 
secreting surfactin divide with probability Pd . The 
newborn cells are placed in neighboring tiles in the 
direction of minimal mechanical resistance. In early 
stages, this may be taken to be the shortest dis¬ 
tance to the biofilm-air interface. 


(dj PPb phase. Accumulation of EPS may prompt 
the creation of an EPS phase according to Eq. ( 12 ). 


(e) Growth tensor and elastic parameters. As EPS 
production increases, the biofilm acquires consis¬ 
tency. We have used in our simulations average 
measured values of the Young modulus. Informa¬ 
tion on the spatial distribution of different types of 
cells allows to spatially modulate that value. Com¬ 
puting the out-of-plane deformations in the next 
stage requires the previous derivation of a growth 
tensor from the biomass growth process. This is 
done as indicated in subsection IIV Cl 


3. Mechanical processes. 

(a) Water absorption and water phase. EPS pro¬ 
duction changes the osmotic pressure and prompts 
water migration from the substratum towards the 
biofilm. Water absorption by the biomass may be 
accounted for by means of formula (27). Dead 


zones near the surface will easily fill with water, 
which may expand along the wrinkles between the 
biomass and the substratum. A water phase is cre¬ 
ated. As the volume fraction of water increases and 
water channels develop, a description of liquid flow 
in the system might be necessary. 


(b) Growth tensor and elastic parameters. Infor¬ 
mation on the spatial distribution of absorbed wa¬ 
ter allows to spatially modulate the average Young 
modulus, decreasing it in swollen regions. The 
growth tensor is updated to reflect the volume sup¬ 
ply due to water absorption. This is done as indi¬ 
cated in subsection IIV Cl 


(c) Out-of-plane displacements. Vertical displace¬ 
ments due to internal stresses are estimated solving 
Eqs. (21)-(22). The equations are nondimension- 
alized following Eqs. ([23])- (24). Basic explicit dif¬ 


ference schemes provide low cost approximations. 
Alternatively, faster spectral methods may be used 
[40 , 42 . As the biofilm becomes more heteroge¬ 
nous we might need to solve a three dimensional 
elasticity problem or include lower order terms in 
the Von-Karman equations. 


(d) Biofilm deformation. Tiles (ix,iy,iz) at the 
interface biofilm-substratum are shifted vertically 
a number of tiles equal to the integer part of 
£(ix, iy)/a , pushing their neighbors in the three di¬ 
mensional biofilm. This is done shifting the status 
of such tiles in the matrix S(ix,iy,iz)- Relative 
status of neighbors in the vertical direction is pre¬ 
served. 
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4. Iteration. The biofilm evolution is calculated al¬ 
ternating the computation of cellular processes up 
to a time tc (estimated from the doubling time), 
with steps of the water absorption processes in a 
time scale tw ~ tc (estimated from biomass pro¬ 
duction) and steps of the deformation processes in a 
time scale tjj (which requires precise measurements 
of the substratum parameters). 

Summarizing, the mechanical processes change the ge¬ 
ometry of the biofilm. The equations for the differ¬ 
ent concentrations must be solved in the new geometry, 
which affects the cellular processes altering the probabili¬ 
ties for the different behaviors. The cellular processes al¬ 
ter the elastic parameters and the pressure in the biofilm, 
contributing residual stresses for its deformation due to 
growth, death and fluid migration. 


VI. SIMULATION RESULTS 


Series of simulations were performed to illustrate the 
behavior of the model and its limitations, choosing sin¬ 
gle mechanisms or combinations of a few of them. The 
results show wrinkled patterns and cell distributions in 
qualitative agreement with recent experimental observa¬ 
tions. They provide insight on the influence of parame¬ 
ters, cellular activities and mechanical processes on the 
biofilm shape and structure. 

Basic growth and spread mechanisms do not produce 
wrinkled shapes. Even if we activate water absorption 
processes, the biofilm will spread faster, but no wrinkles 
appear. When elastic deformation mechanisms are incor¬ 
porated, wrinkles begin to form. We revisit the simula¬ 
tion described in Figure [7] for the same parameter values 
and a similar initial biofilm seed containing five peaks, 
see Figure [llja). During the third step of the growth, 
spread and differentiation processes, EPS producers ap¬ 
pear. At the fourth step, we consider that enough matrix 
has been produced to regard the biofilm as an elastic film 
with Young modulus E. During the 10-th step, dead cells 
appear at the bottom of the initial peaks. Dead areas ex¬ 
pand during the 11-th step. In the 12-th step, growth in 
the central region has slowed down due to depleted nu¬ 
trient levels and increased waste presence. 

Figures [lljb) -(f) illustrate the spatial structure of the 
corresponding growth tensor. Initial peaks diffuse as the 
biofilm grows and become depressions when dead cells ap¬ 
pear. Once growth in the core is depleted, compression is 
higher in the border regions. A randomly perturbed ini¬ 
tial biofilm deforms according to Eqs. (23)-(24) generat¬ 


ing small wrinkles that coarsen as time evolves. Figures 
mb) -(k) represent the out-of-plane deformation of the 
film at selected steps of the growth, spread and differen¬ 
tiation processes. After each step, the average radius of 
the biofilm increases about one cell. The deformation of 
the biofilm is then calculated for a time ^ s. As wrin¬ 
kles develop, the growth processes must be implemented 



FIG. 11. (Color online) (a) Initial biofilm seed, (b), (c), (d), 
(e), (f) Residual stresses e% x averaged after 100 trials of the 
stochastic processes at steps 4, 8 , 10, 11, 12. The average 
compression is —0.132, —0.1425, —0.146, —0.1437, —0.124. 
respectively. £ yy has a similar structure and average. e xy 
takes smaller values and is neglected, (g), (h), (i), (j), (k) 
Vertical displacements. The height mark dx corresponds to 
one cell. Parameter values are 7 = 8 , ^ = 0.5, F = 25 kPa , 
u v — 0.45 and fi v — 0. 


in biofilms with a wavy bottom, contributing additional 
spatial variations to the growth tensor (higher compres¬ 
sion in the valleys, lower compression in the peaks). 

Unless we activate the water absorption mechanism, as 
in Figure [9j a necrotic region will develop in the biofilm 
core. Water absorption facilitates nutrient diffusion and 
waste removal, maintaining the cell normal functions. 

In the computations, we have approximated the resid¬ 
ual stresses £ xx and Sy y by their constant average, mod¬ 
ified by peaks or depressions that represent increased or 
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decreased compression in certain areas. More accurate 
smooth approximations can be automatically produced 
using denoising strategies borrowed from image process¬ 
ing m- As time evolves, singularities may develop due 
to the presence of ‘hanging’ cells in the border. Such cells 
are poorly connected to the rest of the biofilm. This ar¬ 
tifact is solved by discarding such cells when computing 
deformations. For the selected parameter values, biofilms 
are roundish. Therefore, we may alternatively smooth 
the biofilm border using the averaged support obtained 
when averaging trials of the stochastic processes to iden¬ 
tify spatial variations of the growth tensor. Wavy borders 
like the ones observed in Figure [ 6 ] are usually found for 
large values of . This number depends on the nutri¬ 
ent source and t£e bacterial strain, but tends to be small 
( << 1 ) in practice. 

What causes the successive wrinkle branching and the 
radial branching observed in Figures [lja) and[ 2 ja)? As 
we have noticed, EPS production gives the cellular ag¬ 
gregate a certain cohesion. As shown by Figures [To| and 
El an expanding biofilm is under compression due to cell 
division, EPS production and water absorption. The de¬ 
velopment of wrinkles in Figure [Tl] is limited by size con¬ 
siderations. The radius of the biofilm increases from 60 
to 80 cells, it thickness is about 10 cells and 7 = 8 . Dou¬ 
bling the biofilms maximum radius, and consequently 7 , 
it appears that a simple round compression front expand¬ 
ing at a certain speed may produce wrinkle branching, 
structured differently depending on the front speed and 
the compression magnitude. 


Motivated by previous observations, the residual 
stresses £°(x, t) are assumed to behave like a radial front 
e°(r — vt ) expanding at a certain speed u, r being the dis¬ 
tance to the seed center. This allows to lower the compu¬ 
tational cost. Figures |2|b) and [ 8 ] are generated deforming 
circular biofilm seeds according to Eqs. (23)-(|24|) for such 
residual stresses. The evolution starts from a configura¬ 
tion with small random vertical displacements. Figure [ 8 ] 
takes e°(r) = —0.1 constant in the central region. Then, 
it decreases sharply. The front advances one tile every 
14/r s. Figure 12 shows the time evolution of the wrin¬ 
kled area. Small wrinkles similar to those in Figure |TT] 
form that coarsen as shown in those images. Figure [ 8 ] is a 
three dimensional view of the two dimensional projection 
depicted in Figure |T2^d) . The biofilm contains scattered 
dead spots where the residual stresses decrease by a cer¬ 
tain factor, affecting the way wrinkles nucleate. Once 
the wrinkled area attains a certain extension, the height 
of the wrinkles decreases unless we vary the compression 
magnitude. If we wish to increase the wrinkle branching 
rate, while maintaining or enhancing their height as the 
wrinkled region expands, the magnitude of the compres¬ 
sion has to increase radially. Figure [2^b) was computed 
raising the speed to one tile every 1.4 ~fr s and increasing 
the compression magnitude by r/80 as the radius of the 
compressed region grows. When the compression front 
expands too slowly, rings may form around the wrinkled 
area. If it expands too fast, we may see different types 


( a ) f ( b ) £ 



(c) € (d) £ 



FIG. 12. (Color online) Snapshots of the formation of wrin¬ 
kles for £0 = 0 . 1 . (a), (b), (c), (d) images taken at times 

280 — , 560 — , 840 — , 1120 — . An initial random state coarsens 
to produce a connected wrinkle network that opens up into 
radial branches as the compressed region spreads. The biofilm 
has Poisson ratio v — 0.5 and Young modulus E = 25 kPa. 
The Poisson ratio and rubbery modulus of the substratum are 
is v = 0.45 and /i v — 0. T and r are defined in Eq. (25). The 
ratio of the in-plane spatial scale to the out-of-plane spatial 
scale is set to 7 = 16. The equations are nondimensionalized 
so that the dimensionless biofilm thickness becomes h — 1 . 
The time step is dt = ^^dx 2 . The spatial step is dx = 0.1 h. 


of geometric shapes. In between, the network of wrinkles 
opens up radially as it expands. The outer branches are 
connected to the inner network and interact with it. 

So far, we have discussed branching of wrinkles. What 
produces the wrinkled corona in Figure []Ja)? A sud¬ 
den stop of an advancing compression front may arrange 
branching wrinkles in a corona type structure temporar¬ 
ily. The whole network changes its structure later. A uni¬ 
formly compressed film whose elastic modulus decreases 
in an outer corona develops stable radial wrinkles, as in 
Figure [ljb). The Young modulus E has been decreased 
by a factor 0.5 in the outer corona. Radial residual com¬ 
pression of constant magnitude — 0.1 is applied every¬ 
where, for the same parameter values as in Figures [ 8 ] 
and [ 12 ] Now, there is no advancing front. The biofilm 
boundary remains fixed and the residual stresses are con¬ 
stant everywhere. Radial decrease of the Young modulus 
might be due to larger water absorption rates in the outer 
regions, presumably softer, due to the larger presence of 
normal cells. 

Too sharp a decrease in the elastic modulus, for the se¬ 
lected parameters, combined with smaller values of Pois¬ 
son’s ratio, may result in fast coarsening and disappear¬ 
ance of the central wrinkles, as in Figure [l3ja). The 
compression factor is kept constant in a central core but 
increases radially in an outer corona. The inner network 
of wrinkles vanishes with time, but a corona of increas- 
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ing radial wrinkles is formed. Both the Young modulus 
E and the Poisson ratio v affect the time scales. Wrinkles 
coarsen faster in regions with higher Young modulus or 
lower Poisson ratio. Central wrinkles, however, are an¬ 
chored by the presence of a couple of dead spots in Fig¬ 
ure Elb). Decreasing the Poisson ratio enhances height. 
Lowering the Poisson’s ratio to 0.3, the outer wrinkles 
not only maintain their heights, but also may increase it 
for increasing radial compression. Unlike Figure [TJb) or 
Figure [ 8 j a radially increasing compression factor so pro¬ 
duces an inner network of wrinkles that opens up form¬ 
ing branches whose height can be kept uniform over the 
biofilm, see Figure pB^c). Notice also that 7 is decreased 
by a half. Radial increase of the compression might be 
due to larger growth rates in the outer regions caused by 
increasing availability of nutrients. 


VII. CONCLUSIONS 

Replicating cell populations create three dimensional 
organisms of diverse shapes. Elucidating how those dif¬ 
ferent geometries arise is an intriguing question that has 
motivated many theories. Small cellular systems, such as 
bacterial biofilms, seem to develop from the interplay be¬ 
tween cellular and mechanical processes. Studying the in¬ 
teraction of those two mechanisms from an experimental 
point of view requires concurrent measurement of both 
processes, posing a major challenge. Developing mathe¬ 
matical and computational frameworks able to incorpo¬ 
rate the increasing amount of experimental observations 
in their governing rules may provide insight on the feed¬ 
back between microscopic cell behavior and macroscopic 
continuum processes. We have proposed a hybrid model 
to describe the growth dynamics of a small cluster of 
biofilm on an air-agar interface. We are able to trans¬ 
fer information from a stochastic description of cellular 
processes into a continuum model for the deformation 
of the film, and use these deformations as a feedback 
for the cell activities. We have shown that mechanical 
processes change the geometry of the biofilm. The rele¬ 
vant concentrations must be computed in a new geome¬ 
try. This fact influences the cellular processes, modifying 
the probabilities for the cell behaviors. In turn, the cel¬ 
lular processes modulate the elastic parameters and the 
pressure in the biofilm, and generate residual stresses due 
to death, growth, and fluid migration. This microscopic 
feedback determines the subsequent mechanical processes 
and so on. We have observed that the properties of the 
substratum together with the mechanical properties of 
the biofilm and residual stresses due to death, growth 
and swelling seem to control wrinkled biofilm shapes. 

Our in silico biofilms agree qualitatively with some ex¬ 
perimental observations, in the sense that their shape 
and the spatial distribution of the different types of dif¬ 
ferentiated cells is similar. Biomass, air, water and agar 
distribute through the tiles of a computational grid. We 
confer tiles occupied by alive cells the ability to change 



(b) 

1.8 dx 


0 


-1.8 dx 




FIG. 13. (Color online) Wrinkled structures obtained vary¬ 
ing the residual stresses and the elastic moduli for a smaller 
film with lower Poisson’s ratio, (a) Transition from a harder 
central region with larger Young modulus to a softer outer 
corona: radial branches in the corona are separated from the 
core by rings, (b) Uniform residual compression in the cen¬ 
tral region with localized sinks due to dead spots plus radially 
increasing compression in the outer corona: small wrinkles an¬ 
chored by the dead spots in the center merge with larger ra¬ 
dial wrinkles in the outer corona, (c) Residual compression in¬ 
creasing slowly with distance to the center: labyrinths formed 
in the center become radial branches in the outer corona. Pa¬ 
rameter values: Same as previous figures except v — 0.3 and 
7 = 8 . dx is the spatial step. 


their status according to probability laws informed by a 
cascade of concentration fields: nutrients, waste, ComX 
and surfactin. Local inhomogeneities trigger cell death at 
the bottom of biofilm peaks. Nutrient depletion deacti¬ 
vates cells at the top. When a threshold level of ComX is 
reached, surfactin producers appear. When a threshold 
level of surfactin is achieved, EPS producers proliferate. 
Once EPS release activates, we consider the biofilm an 
elastic film, whose elastic parameters are modulated by 
its microstructure. We also launch a mechanism for wa¬ 
ter absorption triggered by variations in osmotic pressure 
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caused by EPS production. Weighting at each biofilm lo¬ 
cation the average values of elastic moduli taken from 
experiments according to the status of the neighboring 
tiles, we may produce spatially varying parameters tak¬ 
ing into account the microstructure. The biofilm weak¬ 
ens due to the presence dead cells. It also softens in 
swollen regions, due to water absorption. Biofilm expan¬ 
sion is quantified by means of a growth tensor computed 
from the stochastic growth, death and water absorption 
processes. In this way, we estimate the residual stresses 
caused by these mechanisms. The compression residual 
stresses and the spatially varying moduli used in the nu¬ 
merical tests presented here are motivated by those com¬ 
putations. We study biofilm deformation in response to 
those stresses by means of a Foppl-Von Karman approx¬ 
imation, obtaining intricate wrinkled cores that split in 
radial branches. Water erodes the regions were cells die, 
expanding along the wrinkles. This improves transport 
of waste and nutrients, hindering cell death and favoring 
growth. 

Our simulations show that wrinkling seems to be as¬ 
sociated with stiffness fluctuations. Variations in inner 
residual stresses and elastic constants due to growth, 
swelling and death seem to govern wrinkle nucleation and 
branching. Compression stresses expanding radially due 
to swelling and growth produce wrinkled cores that split 
in radial wrinkles. The presence of dead regions alters 
the way the core wrinkles are nucleated. It also favors 
wrinkle formation and persistence around dead areas as 
the compression rate is reduced or time increases. Suc¬ 
cessive wrinkle branching may occur depending on the 
expansion velocity. Radially graded compression stresses 
may enhance the outer radial wrinkles. Spatially depen¬ 
dent elastic moduli, harder at the center and softer in the 


outer corona, also enhance the outer radial wrinkles. 

Quantitative comparison with experiments should re¬ 
quire careful calibration of several parameters, in par¬ 
ticular, regarding substratum properties and heuristic 
probability laws. Substratum parameters such as its 
Poisson ratio, thickness, viscosity and rubbery modulus 
should be accurately measured since they affect the time 
scales for the biofilm dynamics. Uptake rates, produc¬ 
tion rates, saturation constants and diffusivities of the 
involved chemical compounds should be determined too, 
since they influence the extend of growth, death and dif¬ 
ferentiation processes. As the biofilm thickens and its 
spatial heterogeneity increases, the Foppl-Von Karman 
approximation should likely be replaced by a fully three 
dimensional elasticity model. The spread mechanism and 
the dynamics of water in the system should be revised 
too. 
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